home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / DBio / DSeqFile.cpp < prev    next >
Text File  |  1995-12-17  |  11KB  |  455 lines
  1. // DSeqFile.cp 
  2. // d.g.gilbert
  3.     
  4.  
  5. #include <ncbi.h>
  6. #include <dgg.h>
  7. #include <stdio.h>
  8. #include <DFile.h>
  9. #include "DSeqFile.h"
  10. #include "DSeqList.h"
  11. #include "ureadseq.h"
  12.  
  13.  
  14.  
  15. Global char* gOutputName = NULL;
  16.  
  17. Boolean DSeqFile::gWriteMasks = true;
  18. static char* kMaskName = "#Mask";
  19.  
  20. #if NEED_LATER
  21. // must be in scope of a DSeqFile:: method
  22. Local const  struct formatTable {
  23.   char  *name;
  24.   short num;
  25.   } formname[] = {
  26.     {"ig",  kIG},
  27.     {"stanford", kIG},
  28.     {"genbank", kGenBank},
  29.     {"gb", kGenBank},
  30.     {"nbrf", kNBRF},
  31.     {"embl", kEMBL},
  32.     {"gcg", kGCG},
  33.     {"uwgcg", kGCG},
  34.     {"dnastrider", kStrider},
  35.     {"strider", kStrider},
  36.     {"fitch", kFitch},
  37.     {"pearson", kPearson},
  38.     {"fasta", kPearson},
  39.     {"zuker", kZuker},
  40.     {"olsen", kOlsen},
  41.     {"phylip", kPhylip},
  42.     {"phylip3.2", kPhylip2},
  43.     {"phylip3.3", kPhylip3},
  44.     {"phylip3.4", kPhylip4},
  45.     {"phylip-interleaved", kPhylip4},
  46.     {"phylip-sequential", kPhylip2},
  47.     {"plain", kPlain},
  48.     {"raw", kPlain},
  49.     {"pir", kPIR},
  50.     {"codata", kPIR},
  51.     {"asn.1", kASN1},
  52.     {"msf", kMSF},
  53.     {"paup", kPAUP},
  54.     {"nexus", kPAUP},
  55.     {"pretty", kPretty},
  56.     {0,0}
  57.   };
  58. #endif
  59.  
  60. static char *gSeqFileFormats[] = { //DSeqFile::kMaxFormat+1
  61.         "IG|Stanford",
  62.         "GenBank|GB",
  63.         "NBRF",
  64.         "EMBL",
  65.         "GCG",
  66.         "DNAStrider",
  67.         "Fitch",
  68.         "Pearson|Fasta",
  69.         "Zuker [in only]",
  70.         "Olsen [in only]",
  71.         "Phylip3.2|Phylip2",
  72.         "Phylip|Phylip4",
  73.         "Plain|Raw",
  74.         "PIR|CODATA",
  75.         "MSF",
  76.         "ASN.1",
  77.         "PAUP|NEXUS",
  78.         "Pretty [out only]",
  79.         0};
  80.  
  81. const char* DSeqFile::FormatName(short index)
  82. { 
  83.     if (index<1 || index>18) return "";
  84.     else return gSeqFileFormats[index-1];
  85. }
  86.  
  87. short  DSeqFile::FormatFromName(const char* name)
  88. {
  89.     if (!name) return kUnknown;
  90.     long   i, len1, namelen= StrLen(name);
  91.     char      *form1, *form2;
  92.     if (namelen<2) return kUnknown;
  93.     for (i=0; i<kMaxFormat; i++) {
  94.         form1= (char*) gSeqFileFormats[i];
  95.         form2= StrChr(form1,'|'); 
  96.         if (form2) { len1= form2 - form1; form2++; }
  97.         else len1= StrLen( form1);
  98. #if 1
  99.         if (StrNICmp(name, form1,len1) == 0) return i+1;
  100.         else if (form2 && StrICmp(name, form2) == 0) return i+1;
  101. #else
  102.             // problem w/ Phylip matching Phylip3.2 not Phylip
  103.         if (StrNICmp(name, form1, namelen) == 0) return i+1;
  104.         else if (form2 && StrNICmp(name, form2, namelen) == 0) return i+1;
  105. #endif
  106.         }
  107.     return kUnknown;
  108. }
  109.  
  110.  
  111. char* DSeqFile::FixSeqID( char* s, short leftmargin)
  112. {
  113. #if 1
  114.     register char * cp, * ep;
  115.     for (cp= s; *cp && isspace(*cp); cp++) ;
  116.     for (ep= cp; *ep && !isspace(*ep); ep++) ;
  117.     if (ep>cp && ep[-1] == ',') ep--; // fix for fasta "name, #bases, # checksum"
  118.     if (*ep) *ep= 0;
  119.     return cp;
  120. #else
  121.     //s= StrDup(s);
  122.   short len = StrLen(s);
  123.     while (s[--len] == ' ' &&  len>0) ;
  124.     s[++len]= 0;
  125.     short i = Min( len-1, leftmargin);
  126.     while (s[i] == ' ' && i<len) i++;
  127.     if (i>0) s += i; //MemMove(s, s+i, len-i+1); //s= s+i;
  128.     //if (StrLen(s) > kNameWidth) s[kNameWidth]= 0;
  129.     return s;
  130. #endif
  131. }
  132.  
  133. short DSeqFile::SeqFileFormatWrapper( DFile* aFile) 
  134. {
  135.     long    skiplines = 0;
  136.     short    aFormat, err = 0;
  137.  
  138.     aFile->Open("r");  
  139.     aFormat= seqFileFormatFp( aFile->fFile, &skiplines, &err);
  140.     aFile->Close();  
  141.     return    aFormat;
  142. }
  143.  
  144.  
  145. short DSeqFile::ReadSeqFile( DFile* aFile, DSeqList* aSeqList) 
  146. {
  147.     char    seqid[256], *seqidptr, *bases, *seqlist;
  148.     long    skiplines = 0, skipi, seqlen;
  149.     short aFormat, atseq = 0, whichSeq = 0, nseq = 0, err = 0;
  150.     char  *defname;
  151.     Boolean notdone = true, needrewind = false;
  152.     DSequence* aSeq;
  153.     
  154.     aFile->OpenFile("r");    
  155.     aFormat = seqFileFormatFp( aFile->fFile, &skiplines, &err);
  156.     aFile->OpenFile(); // rewind
  157.     defname= FixSeqID( gDefSeqName, 0);
  158.     
  159.     switch (aFormat) {
  160.         case kUnknown:
  161.         case kNoformat:
  162.             aFile->CloseFile();
  163.             return aFormat;
  164.          
  165.             // all multi-pass (interleaved) formats go here !        
  166.         case kOlsen:
  167.         case kMSF:
  168.         case kPAUP:
  169.         case kPhylip2:
  170.         case kPhylip4:
  171.             needrewind = true;
  172.             aFile->CloseFile();
  173.             seqlist = listSeqs( aFile->GetName(), skiplines, aFormat, &nseq, &err);
  174.             MemFree(seqlist); 
  175.             break;
  176.         }
  177.         
  178.     skipi= skiplines;
  179.     while (notdone) {
  180.             // this one-pass thru gFileRef may well be buggy for current readSeq
  181.             // need to rewind file & count whichSeq++ for each sequence...?
  182.         if (needrewind) {
  183.             whichSeq++; 
  184.             atseq= 0; 
  185.             aFile->OpenFile("r");
  186.             skipi= skiplines;
  187.             }
  188.         else {
  189.             whichSeq = 1; 
  190.             atseq= 0;
  191.             }
  192.         seqlen = 0; seqidptr= seqid; *seqidptr= 0;
  193.         bases = readSeqFp( whichSeq, aFile->fFile, skipi, aFormat, &seqlen, 
  194.                                                         &atseq, &err, seqidptr);
  195.         skipi= 0; 
  196.         
  197.                 // install bases into new TSequence, and stick TSequence into list
  198.         if (bases) {
  199.             if (seqlen > 0) {
  200.                 char *name, * cp;
  201.                 Boolean ismask= false;
  202.                 
  203.                 if (*seqidptr != 0) {
  204.                     if ( (cp= StrStr(seqidptr, kMaskName)) != NULL) {
  205.                         ismask= true;
  206.                         *cp= 0;
  207.                         }
  208.                     name= FixSeqID( seqidptr, 0);
  209.                     }
  210.                 else 
  211.                     name= defname;
  212.                     
  213.                 if (ismask) {
  214.                     aSeq= aSeqList->FindNamedSeq( name);
  215.                     if (aSeq) {
  216.                         // convert mask...
  217.                         register char *mp;
  218.                         for (mp= bases; *mp; mp++) { 
  219.                             *mp -= '?'; // to zero relative
  220.                             *mp |= 0x80; // add high bit
  221.                             }
  222.                         aSeq->SetMasks( bases);
  223.                         }
  224.                     }
  225.                 else {
  226.                     aSeq = MakeSequence( name, bases, NULL, 0); //gFileRef->moddate
  227.                     aSeqList->InsertLast(aSeq);
  228.                     }
  229.                 }
  230.             else notdone= false;
  231.             MemFree(bases); // readSeq made this ptr w/ malloc call... MakeSeq copies this data... 
  232.             }
  233.         else notdone= false; // nodata -- if format=0, no data is read
  234.             
  235.         if (needrewind) notdone &= (whichSeq < nseq); // stop now
  236.         else notdone &= !aFile->Eof();  
  237.         }
  238.         
  239.     aFile->CloseFile();
  240.     return aFormat;
  241. }
  242.  
  243.  
  244.  
  245.  
  246. void    DSeqFile::WriteSeqHeader( DFile* aFile,  
  247.                                 short format, short nseqs, long nbases, char* outputName,
  248.                                 Boolean& needSameSize, Boolean& isInterleaved)
  249. {
  250.     const char *kASN1headline = "Bioseq-set ::= {\nseq-set {\n";
  251.     char    line[256], name[64];
  252.     
  253.     // this is a ureadseq.h call
  254.     // need to do somewhere once, then let user set prefs...
  255.  
  256. #define gPrettyInit1(p) { \
  257.   p.isactive=false;\
  258.   p.baseonlynum=true;\
  259.   p.numline= p.atseq= 0;\
  260.   p.numright= p.numtop= true; \
  261.   p.numleft= p.numbot= false;\
  262.   p.nameright= true; \
  263.   p.nameleft= p.nametop= false;\
  264.   p.noleaves= p.domatch= p.degap= false;\
  265.   p.matchchar='.';\
  266.   p.gapchar='-';\
  267.   p.namewidth=8;\
  268.   p.numwidth=5;\
  269.   p.interline=1;\
  270.   p.spacer=10;\
  271.   p.seqwidth=50;\
  272.   p.tab=0; }
  273.  
  274.     gPrettyInit1(gPretty);
  275.     //gPretty.interline= 2; // not active, along w/ numtop,numbot,nametop...
  276.     // need to work into this class
  277.     
  278.     
  279.         
  280.     needSameSize= false;
  281.     isInterleaved= false;
  282.     aFile->Open("a"); 
  283.  
  284.     if (outputName) StrNCpy( name, outputName, 63);
  285.     else StrCpy( name, DFileManager::kUntitled);
  286.     
  287.     switch (format) {
  288.         case kPIR:
  289.             aFile->WriteLine( "\\\\\\\n");
  290.             break;
  291.         case kASN1:
  292.             aFile->WriteLine( (char*)kASN1headline);
  293.             break;
  294.  
  295.         case kMSF:     // this is for top of interleaved file
  296.             {
  297.             isInterleaved= true;
  298.          unsigned long checksum = 0, checkall = 0;
  299.             //checksum= GCGchecksum(seq, seqlen, &checkall);
  300.             sprintf(line,"\n %s  MSF: %d  Type: N  June 01, 1776 12:00  Check: %d ..\n\n",
  301.                                         name, nbases, checkall);
  302.             aFile->WriteLine( line);
  303.             }
  304.             break;
  305.  
  306.         case kPhylip2:     
  307.             isInterleaved= false;
  308.             goto allphyl;
  309.             break;
  310.  
  311.         case kPhylip4:      
  312.             isInterleaved= true;
  313.         allphyl:
  314.             {
  315. //    if (phylvers >= 4) fprintf(foo," %d %d\n", seqout, seqlen);
  316. //    else fprintf(foo," %d %d YF\n", seqout, seqlen);
  317.             if (isInterleaved) sprintf( line, " %d %d F\n", nseqs, nbases);  
  318.             else sprintf( line, " %d %d FI \n", nseqs, nbases);  
  319.             aFile->WriteLine( line);
  320.             needSameSize= true;
  321.             }
  322.             break;
  323.  
  324.         case kPAUP:
  325.             isInterleaved= true;
  326.             needSameSize= true;
  327.             aFile->WriteLine( "#NEXUS\n");
  328.             sprintf( line,"[%s -- data title]\n\n", name);
  329.             aFile->WriteLine( line);
  330.             break;
  331.  
  332.         case kPretty:
  333.             isInterleaved= true;
  334.             gPretty.isactive= true;
  335.             break;
  336.             
  337.         default:
  338.             break;
  339.         }
  340.         
  341.     aFile->CloseFile();
  342. }
  343.  
  344.  
  345.  
  346. void    DSeqFile::DoInterleave( DFile* toFile, DFile* fromFile, 
  347.                                 short outform, short seqtype, short nseqs, long nbases,
  348.                                 short nlines, short noutindex, long* outindex)
  349. {
  350.     short        leaf, iline, iseq;
  351.     char        *cp, line[256];
  352.     
  353.     toFile->Open("a");
  354.     fromFile->Open("r");
  355.     
  356.     for (leaf=0; leaf<nlines; leaf++) {
  357.     
  358.         if (outform == kMSF && leaf == 1) { // do after writing seq names
  359.             toFile->WriteLine( "//\n\n");
  360.             }
  361.             
  362.         if (outform == kPAUP && leaf==1) { // do after writing seq names
  363.             switch (seqtype) {
  364.                 case DSequence::kDNA     : cp= "dna"; break;
  365.                 case DSequence::kRNA     : cp= "rna"; break;
  366.                 case DSequence::kNucleic : cp= "dna"; break;
  367.                 case DSequence::kAmino   : cp= "protein"; break;
  368.                 case DSequence::kOtherSeq: cp= "dna"; break;
  369.                 }
  370.             toFile->WriteLine( "\nbegin data;\n");
  371.             sprintf(line," dimensions ntax=%d nchar=%d;\n", nseqs, nbases); 
  372.             toFile->WriteLine( line);
  373.             sprintf(line," format datatype=%s interleave missing=%c", cp, gPretty.gapchar); // gapchar == indelHard ??? 
  374.             toFile->WriteLine( line);
  375.             //if (gPretty.domatch) { sprintf(line," matchchar=%c", gPretty.matchchar); 
  376.                 //Mwrite( toFile, line); }
  377.             toFile->WriteLine( ";\n  matrix\n");
  378.             }
  379.  
  380.         for (iseq=0; iseq<noutindex; iseq++) {
  381.             fromFile->SetDataMark( outindex[iseq]);
  382.             
  383.             for (iline=0; iline<=leaf; iline++)
  384.                 if (fromFile->ReadLine( line, 256) != 0) *line= 0;
  385.                 
  386.             if (fromFile->Tell() <= outindex[iseq+1])
  387.                 toFile->WriteLine( line, false);  //remember- ReadLine a.k.a fgets retains newline
  388.             }
  389.  
  390.      for (iline=0; iline<gPretty.interline; iline++)
  391.         toFile->WriteLine("\n"); //LineEnd 
  392.             
  393.         }
  394.         
  395.     fromFile->Close();
  396.     toFile->Close();
  397. }
  398.  
  399.  
  400. void    DSeqFile::WriteSeqTrailer( DFile* aFile,
  401.                                 short format, short /*nseqs*/, short /*nlinesperseq*/, long /*nbases*/)
  402. {
  403.     switch (format) {
  404.   
  405.       case kASN1:  
  406.             aFile->Open("a"); 
  407.             aFile->WriteLine("} }\n"); 
  408.             aFile->Close(); 
  409.             break;
  410.     
  411.       case kPAUP:    
  412.             aFile->Open("a"); 
  413.             aFile->WriteLine(";\n  end;\n");
  414.             aFile->Close(); 
  415.             break;    
  416.         }
  417. }
  418.  
  419.  
  420.  
  421. void DSeqFile::WriteSeqWrapper( DFile* aFile,
  422.                                     char *seq, long seqlen, short outform, char *seqid)
  423. { 
  424.     if (!aFile->IsOpen()) aFile->Open("a");
  425.     seqid= StrDup(seqid);
  426.     char* cp= FixSeqID(seqid, 0);
  427.     gLinesPerSeqWritten= writeSeq( aFile->fFile, seq, seqlen, outform, cp);
  428.     MemFree(seqid);
  429.     //aFile->Close(); //?? want left open for appending by caller ? e.g., seqmail
  430. }  
  431.  
  432. void DSeqFile::WriteMaskWrapper( DFile* aFile,
  433.                                     char *mask, long seqlen, short outform, char *seqid)
  434. { 
  435.     char maskid[128];
  436.     if (!aFile->IsOpen()) aFile->Open("a"); 
  437.     seqid= StrDup(seqid);
  438.     char* cp= FixSeqID(seqid, 0);
  439.     StrNCpy( maskid, cp, sizeof(maskid)-5);
  440.     StrNCat( maskid, kMaskName, sizeof(maskid));
  441.     MemFree(seqid);
  442.         //     massage mask into writable form
  443.     mask= StrDup(mask); 
  444.     Boolean hasdata= false;
  445.     for (cp= mask; *cp; cp++) { 
  446.         *cp &= 0x7F; // drop hi bit
  447.         if (*cp != 0) hasdata= true;
  448.         *cp += '?';  // set 0 == '?', mask range is mostly 'A-Za-z'
  449.      }
  450.     if (hasdata)
  451.         gLinesPerSeqWritten= writeSeq( aFile->fFile, mask, seqlen, outform, maskid);
  452.     MemFree(mask);
  453. }  
  454.  
  455.